#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdbool.h>

#define MAX_POLY_LEN 100  // 最大多项式系数长度
#define MAX_STR_LEN 1024  // 最大字符串长度

// ========================= 多项式化简工具（复刻polynomial_simplify） =========================
// 多项式乘法：(a0+a1x+...)(b0+b1x+...) → 结果系数数组
double* poly_multiply(double* a, int a_len, double* b, int b_len, int* out_len) {
    *out_len = a_len + b_len - 1;
    double* res = (double*)calloc(*out_len, sizeof(double));
    for (int i = 0; i < a_len; i++) {
        for (int j = 0; j < b_len; j++) {
            res[i + j] += a[i] * b[j];
        }
    }
    return res;
}

// 多项式化简：将因子矩阵（如[[1,-2],[1,-4]]）展开为单个多项式
double* polynomial_simplify(double** matrix, int rows, int cols, int* out_len) {
    if (rows == 0) {
        *out_len = 1;
        double* res = (double*)malloc(sizeof(double));
        res[0] = 1.0;
        return res;
    }

    // 初始多项式为第一个因子
    double* res = (double*)malloc(cols * sizeof(double));
    memcpy(res, matrix[0], cols * sizeof(double));
    int res_len = cols;

    // 逐个因子相乘
    for (int i = 1; i < rows; i++) {
        double* factor = matrix[i];
        int factor_len = cols;
        int new_len;
        double* new_res = poly_multiply(res, res_len, factor, factor_len, &new_len);

        free(res);
        res = new_res;
        res_len = new_len;
    }

    // 去除末尾零（精度优化）
    while (res_len > 1 && fabs(res[res_len - 1]) < 1e-10) {
        res_len--;
    }

    *out_len = res_len;
    return res;
}

// ========================= 拉格朗日分式结构体（LagrangeFraction） =========================
typedef struct {
    double* numerators;  // 分子因子数组（x - xi 的 xi 集合）
    int num_count;       // 分子因子数量
    double denominator;  // 分母（所有(xj-xi)的乘积）
} LagrangeFraction;

// 创建拉格朗日分式
LagrangeFraction* lagrange_fraction_create(double* nums, int num_cnt, double denom) {
    LagrangeFraction* frac = (LagrangeFraction*)malloc(sizeof(LagrangeFraction));
    frac->numerators = (double*)malloc(num_cnt * sizeof(double));
    memcpy(frac->numerators, nums, num_cnt * sizeof(double));
    frac->num_count = num_cnt;
    frac->denominator = denom;
    return frac;
}

// 释放拉格朗日分式内存
void lagrange_fraction_destroy(LagrangeFraction* frac) {
    if (!frac) return;
    free(frac->numerators);
    free(frac);
}

// ========================= 拉格朗日项结构体（LagrangeItem） =========================
typedef struct {
    double y;                      // 数据点的y值
    LagrangeFraction* fraction;    // 对应的分式
} LagrangeItem;

// 创建拉格朗日项
LagrangeItem* lagrange_item_create(double y, LagrangeFraction* frac) {
    LagrangeItem* item = (LagrangeItem*)malloc(sizeof(LagrangeItem));
    item->y = y;
    item->fraction = frac;
    return item;
}

// 释放拉格朗日项内存
void lagrange_item_destroy(LagrangeItem* item) {
    if (!item) return;
    lagrange_fraction_destroy(item->fraction);
    free(item);
}

// ========================= 拉格朗日多项式结构体（LagrangePolynomial） =========================
typedef struct {
    LagrangeItem** items;          // 拉格朗日项数组
    int item_count;                // 项数量
    double* polynomial;            // 化简后的多项式系数（霍纳形式）
    int poly_len;                  // 多项式系数长度
} LagrangePolynomial;

// 创建拉格朗日多项式（核心：化简多项式）
LagrangePolynomial* lagrange_poly_create(LagrangeItem** items, int item_cnt) {
    LagrangePolynomial* poly = (LagrangePolynomial*)malloc(sizeof(LagrangePolynomial));
    poly->items = items;
    poly->item_count = item_cnt;
    poly->poly_len = item_cnt;
    poly->polynomial = (double*)calloc(poly->poly_len, sizeof(double));

    // 步骤1：逐个项展开并化简
    for (int i = 0; i < item_cnt; i++) {
        LagrangeItem* item = items[i];
        LagrangeFraction* frac = item->fraction;

        // 计算系数：y / 分母
        double coeff = item->y / frac->denominator;

        // 构建因子矩阵：[[1, -xi] ...]
        int factor_rows = frac->num_count;
        int factor_cols = 2;
        double** factor_matrix = (double**)malloc(factor_rows * sizeof(double*));
        for (int j = 0; j < factor_rows; j++) {
            factor_matrix[j] = (double*)malloc(factor_cols * sizeof(double));
            factor_matrix[j][0] = 1.0;
            factor_matrix[j][1] = -frac->numerators[j];
        }

        // 第一个因子乘以系数（优化：避免后续乘法）
        if (factor_rows > 0) {
            factor_matrix[0][0] *= coeff;
            factor_matrix[0][1] *= coeff;
        } else {
            // 无因子（常数项）
            coeff = item->y / frac->denominator;
        }

        // 多项式化简
        int simplify_len;
        double* simplify_poly = polynomial_simplify(factor_matrix, factor_rows, factor_cols, &simplify_len);

        // 合并到总多项式
        for (int j = 0; j < simplify_len; j++) {
            if (j < poly->poly_len) {
                poly->polynomial[j] += simplify_poly[j];
            }
        }

        // 释放临时内存
        free(simplify_poly);
        for (int j = 0; j < factor_rows; j++) {
            free(factor_matrix[j]);
        }
        free(factor_matrix);
    }

    return poly;
}

// 释放拉格朗日多项式内存
void lagrange_poly_destroy(LagrangePolynomial* poly) {
    if (!poly) return;
    for (int i = 0; i < poly->item_count; i++) {
        lagrange_item_destroy(poly->items[i]);
    }
    free(poly->items);
    free(poly->polynomial);
    free(poly);
}

// 霍纳法则求值（核心优化）
double lagrange_poly_evaluate(const LagrangePolynomial* poly, double x) {
    double result = 0.0;
    for (int i = 0; i < poly->poly_len; i++) {
        result = result * x + poly->polynomial[i];
    }
    return result;
}

// 多项式表达式字符串（复刻__str__）
char* lagrange_poly_to_str(const LagrangePolynomial* poly) {
    char* str = (char*)malloc(MAX_STR_LEN * sizeof(char));
    str[0] = '\0';
    bool first = true;

    for (int i = 0; i < poly->item_count; i++) {
        LagrangeItem* item = poly->items[i];
        LagrangeFraction* frac = item->fraction;

        // 构建分子部分：(x-x1)(x-x2)...
        char num_str[MAX_STR_LEN] = "";
        for (int j = 0; j < frac->num_count; j++) {
            char term[32];
            sprintf(term, "(x-%.0f)", frac->numerators[j]);
            strcat(num_str, term);
        }

        // 构建项字符串：(x-x1)(x-x2)*y/(denom)
        char item_str[MAX_STR_LEN];
        sprintf(item_str, "%s*%.1f/(%.0f)", num_str, item->y, frac->denominator);

        // 拼接符号
        if (!first) {
            strcat(str, "+");
        }
        first = false;
        strcat(str, item_str);
    }

    return str;
}

// ========================= 拉格朗日插值核心函数 =========================
// 创建单个拉格朗日分式（对应create_fraction）
LagrangeFraction* create_fraction(double** data, int data_rows, int j) {
    double* numerators = (double*)malloc((data_rows - 1) * sizeof(double));
    int num_cnt = 0;
    double denominator = 1.0;
    double xj = data[j][0];

    for (int i = 0; i < data_rows; i++) {
        if (i == j) continue;
        double xi = data[i][0];
        numerators[num_cnt++] = xi;
        denominator *= (xj - xi);
    }

    return lagrange_fraction_create(numerators, num_cnt, denominator);
}

// 拉格朗日插值主函数（对应interpolate）
LagrangePolynomial* lagrange_interpolate(double** data, int data_rows) {
    // 构建拉格朗日项数组
    LagrangeItem** items = (LagrangeItem**)malloc(data_rows * sizeof(LagrangeItem*));
    for (int i = 0; i < data_rows; i++) {
        LagrangeFraction* frac = create_fraction(data, data_rows, i);
        items[i] = lagrange_item_create(data[i][1], frac);
    }

    // 创建并化简多项式
    LagrangePolynomial* poly = lagrange_poly_create(items, data_rows);
    return poly;
}

// ========================= 测试代码 =========================
int main() {
    // 测试数据：[[0,1],[2,5],[4,17]]
    int data_rows = 3;
    double** data = (double**)malloc(data_rows * sizeof(double*));
    data[0] = (double*)malloc(2 * sizeof(double));
    data[1] = (double*)malloc(2 * sizeof(double));
    data[2] = (double*)malloc(2 * sizeof(double));
    data[0][0] = 0.0; data[0][1] = 1.0;
    data[1][0] = 2.0; data[1][1] = 5.0;
    data[2][0] = 4.0; data[2][1] = 17.0;

    // 拉格朗日插值计算
    LagrangePolynomial* poly = lagrange_interpolate(data, data_rows);
    if (!poly) {
        fprintf(stderr, "插值失败\n");
        return -1;
    }

    // 测试求值（霍纳法则）
    printf("poly(0) = %.1f, poly(2) = %.1f, poly(4) = %.1f\n",
           lagrange_poly_evaluate(poly, 0.0),
           lagrange_poly_evaluate(poly, 2.0),
           lagrange_poly_evaluate(poly, 4.0));

    printf("poly(1) = %.1f, poly(3) = %.1f, poly(5) = %.1f\n",
           lagrange_poly_evaluate(poly, 1.0),
           lagrange_poly_evaluate(poly, 3.0),
           lagrange_poly_evaluate(poly, 5.0));

    // 打印原始多项式表达式
    char* poly_str = lagrange_poly_to_str(poly);
    printf("原始拉格朗日表达式：%s\n", poly_str);

    // 打印化简后的多项式系数
    printf("化简后多项式系数（霍纳形式）：[");
    for (int i = 0; i < poly->poly_len; i++) {
        printf("%.3f", poly->polynomial[i]);
        if (i != poly->poly_len - 1) {
            printf(", ");
        }
    }
    printf("]\n");

    // 释放内存
    free(poly_str);
    lagrange_poly_destroy(poly);
    for (int i = 0; i < data_rows; i++) {
        free(data[i]);
    }
    free(data);

    return 0;
}